import numpy as np
import matplotlib.pyplot as plt
#绘制三维图像
import mpl_toolkits.mplot3d as p3d

'''
时滞混沌系统Mackey-Glass，Euler(欧拉)公式进行求解,其他方法也可。。。
'''
def Mackey_Glass(x,h,a,b,c,tao,T):
  #这里t是迭代周期，公式x(t)中的t是t*h
  for t in range(tao,tao+T):
    x[t+1]=x[t]+h*(-b*x[t]+a*x[t-tao]/(1+x[t-tao]**c))
  return x

def main():
  #设定参数
  a=0.2
  b=0.1
  c=10
  tao=17
  #迭代次数
  T=1177
  #初值是用连续函数定义的
  #迭代步长，应满足tao%h=0,否则存在误差
  h=1
  #tao转成离散序列索引
  tao=int(tao/h)
  #前tao+1个数是[-tao,0]的系统初值，后T个数是迭代T次
  #dtype默认是int32，要指定dtype
  x=np.arange(0,T+tao+1,dtype=np.float64)
  x1=np.arange(0,T+tao+1,dtype=np.float64)
  x2=np.arange(0,T+tao+1,dtype=np.float64)
  #[-tao,0]之间的初始条件
  for i in range(tao+1):
    x[i]=np.cos(i)
    #对原初始函数加10的负8次方扰动，对比结果
    x1[i]=np.cos(i)+10**(-8)
    x2[i]=np.cos(i)+10**(-8)
    # x[i]=1.2
    # #对原初始函数加10的负8次方扰动，对比结果
    # x1[i]=1.2+10**(-8)
    # x2[i]=1.2+10**(-8)
  #产生时滞混沌序列
  x=Mackey_Glass(x,h,a,b,c,tao,T)
  x1=Mackey_Glass(x1,h,a,b,c,tao,T)

  #结果展示
  # plt.rcParams['font.sans-serif'] = ['SimHei']  # 中文乱码
  # plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
  # plt.subplot(221)
  # #舍去序列中前300个暂态点，这里采用t的值而不是索引
  # t=np.arange(300,T+tao+1)
  # plt.scatter(t,x[300:],s=1)
  # plt.title('$x_0(t)=cos(t),t\\in[-\\tau,0]$')
  # #子图2
  # plt.subplot(222)
  # plt.scatter(t,x1[300:],s=1)
  # plt.title('$x_0(t)=cos(t)+10^{-8},t\\in[-\\tau,0]$')
  # plt.subplot(223)
  # plt.scatter(t,x[300:]-x1[300:],s=1)
  # plt.title('$x_0(t)$与$x_1(t)$的差值')
  # plt.subplot(224)
  # plt.plot(x[300:T],x[300+tao:tao+T])
  # plt.title('x(t)与x(t+tao)的相图')
  #
  # # #设置子图默认的间距
  # plt.tight_layout()
  # #显示图像
  # plt.show()


  plt.figure(figsize=(16, 4))
  plt.plot(np.arange(801,T+tao+1), x[801:], 'r-', label='Test Target')

  plt.xlabel('X')
  plt.ylabel('Y')
  # 设置 y 轴范围
  # plt.ylim(0.4, 1.5)  # y 轴的范围设定为 -2 到 2
  plt.legend()
  plt.show()
if __name__=='__main__':
  main()

